home *** CD-ROM | disk | FTP | other *** search
/ El Mac 9 / El Mac 9.iso / Shareware / Applications / MathPad 2.4 / XFuns / XFun kit / linfit src / linfit68K.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-01-10  |  1.7 KB  |  73 lines  |  [TEXT/CWIE]
  1. /*
  2.  linfit(XYpoints)
  3.  Fit data to a straight line using linear regression formula.
  4. */
  5.  
  6. /* This file uses the "no globals" version of callbacks.
  7.    (PC relative strings). */
  8.    
  9. #include "callback.h"        /* use "no global" version of callbacks */
  10.  
  11. /* can't use math library sqrt() without A4 globals so implement a sqrt() here
  12.    via  FixMath toolbox call. Works only for  .5  < x < 1e6 */
  13. double sqrt(double x);
  14. double sqrt(double x)
  15. {
  16.    return 1/Frac2X(FracSqrt(X2Frac(1/x)));
  17. }
  18.  
  19. static short linfit(double *retval,funptr callback)
  20. {
  21.    EXPR arr;
  22.    double x,y,*iptr,*jptr,intercept;
  23.    double sumx,sumxx,sumy,sumyy,sumxy,Sxx,Sxy;
  24.    short isarray;
  25.    long n,count;
  26.    
  27.    sumx=0; sumxx=0; sumy=0; sumyy=0; sumxy=0;
  28.    n = 0;
  29.    MakeParmExpr(0,&arr,callback);
  30.    ProbeExpr(arr,&x,&isarray,&count,callback);
  31.    if(!isarray || count<2)
  32.    {
  33.     ErrMsg(" linfit() expects array of {x,y}",0L,callback);
  34.     FreeExpr(arr,callback);
  35.     return(FALSE);
  36.    }
  37.    AddIndex(&arr,&iptr,callback);    // arr[i] is {x,y}
  38.    AddIndex(&arr,&jptr,callback);    // j picks x or y
  39.    while(count--)
  40.    {
  41.     *jptr = 1;
  42.     if(EvalExpr(arr,&x,callback) && x==x)
  43.     {
  44.      *jptr = 2;
  45.      if(EvalExpr(arr,&y,callback) && y==y)
  46.      {
  47.       sumx += x;
  48.       sumxx += x*x;
  49.       sumy += y;
  50.       sumyy += y*y;
  51.       sumxy += x*y;
  52.       n++;
  53.      }
  54.     }
  55.     *iptr += 1;
  56.    }
  57.    Sxx=n*sumxx-sumx*sumx;
  58.    Sxy=n*sumxy-sumx*sumy;
  59.  
  60.    SetVarVal("slope",Sxy/Sxx,callback);
  61.    intercept=(sumxx*sumy-sumx*sumxy)/(n*sumxx-sumx*sumx);
  62.    SetVarVal("intercept",intercept,callback);
  63.    *retval = Sxy/sqrt(Sxx*(n*sumyy-sumy*sumy));        // return correlation coeff
  64.    
  65.    FreeExpr(arr,callback);
  66.    return(TRUE);
  67. }
  68.  
  69. void main(funptr callback)
  70. {
  71.     AddXfun("linfit","XYpoints",&linfit,NULL,callback);
  72. }
  73.